<html>
 <head>
  <meta charset="utf-8"/>
  <meta content="width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no" name="viewport"/>
  <title>
   Editor: MCMC案例学习  | 数螺 | NAUT IDEA
  </title>
  <link href="http://cdn.bootcss.com/bootstrap/3.3.6/css/bootstrap-theme.min.css" rel="stylesheet"/>
  <link href="http://cdn.bootcss.com/bootstrap/3.3.6/css/bootstrap.min.css" rel="stylesheet"/>
  <style type="text/css">
   #xmain img {
                  max-width: 100%;
                  display: block;
                  margin-top: 10px;
                  margin-bottom: 10px;
                }

                #xmain p {
                    line-height:150%;
                    font-size: 16px;
                    margin-top: 20px;
                }

                #xmain h2 {
                    font-size: 24px;
                }

                #xmain h3 {
                    font-size: 20px;
                }

                #xmain h4 {
                    font-size: 18px;
                }


                .header {
	           background-color: #0099ff;
	           color: #ffffff;
	           margin-bottom: 20px;
	        }

	        .header p {
                  margin: 0px;
                  padding: 10px 0;
                  display: inline-block;  
                  vertical-align: middle;
                  font-size: 16px;
               }

               .header a {
                 color: white;
               }

              .header img {
                 height: 25px;
              }
  </style>
  <script src="http://cdn.bootcss.com/jquery/3.0.0/jquery.min.js">
  </script>
  <script src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript">
   MathJax.Hub.Config({
          tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]},
          TeX: {equationNumbers: {autoNumber: "AMS"}}
        });
  </script>
  <script src="http://nautstatic-10007657.file.myqcloud.com/static/css/readability.min.js" type="text/javascript">
  </script>
  <script type="text/javascript">
   $(document).ready(function() {
                 var loc = document.location;
                 var uri = {
                  spec: "http://cos.name/2012/07/mcmc-case-study/",
                  host: "http://cos.name",
                  prePath: "http://cos.name",
                  scheme: "http",
                  pathBase: "http://cos.name/"
                 };
    
                 var documentClone = document.cloneNode(true);
                 var article = new Readability(uri, documentClone).parse();
     
                 document.getElementById("xmain").innerHTML = article.content;
                });
  </script>
  <!-- 1466441765: Accept with keywords: (title(0.2):中国,案例,服务平台,门户网站,统计学, topn(0.7):概率,方差,标准,R语言,参数估计,标准差,长度,后验,结果,模型,均值,参数,相关性,计算,数据,正态分布,函数,阶矩,样本,问题,案例,算法,先验,线性,对数,利用,方法,估计值,分布,广义).-->
 </head>
 <body class="single single-post postid-6007 single-format-standard sidebar" onload="">
  <div class="header">
   <div class="container">
    <div class="row">
     <div class="col-xs-6 col-sm-6 text-left">
      <a href="/databee">
       <img src="http://nautidea-10007657.cos.myqcloud.com/logo_white.png"/>
      </a>
      <a href="/databee">
       <p>
        数螺
       </p>
      </a>
     </div>
     <div class="hidden-xs col-sm-6 text-right">
      <p>
       致力于数据科学的推广和知识传播
      </p>
     </div>
    </div>
   </div>
  </div>
  <div class="container text-center">
   <h1>
    Editor: MCMC案例学习
   </h1>
  </div>
  <div class="container" id="xmain">
   <div class="hfeed site" id="page">
    <header class="site-header" id="masthead" role="banner">
     <div id="cos-logo">
      <a href="http://cos.name/">
       <img src="http://cos.name/wp-content/themes/COS-kermesinus/images/headers/cos-logo.png"/>
      </a>
     </div>
     <div class="navbar" id="navbar">
      <nav class="navigation main-navigation" id="site-navigation" role="navigation">
       <h3 class="menu-toggle">
        菜单
       </h3>
       <a class="screen-reader-text skip-link" href="http://cos.name/2012/07/mcmc-case-study/#content" title="跳至内容">
        跳至内容
       </a>
       <div class="menu-%e6%88%91%e7%9a%84%e8%8f%9c%e5%8d%95-container">
        <ul class="nav-menu" id="menu-%e6%88%91%e7%9a%84%e8%8f%9c%e5%8d%95">
         <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-home menu-item-4746" id="menu-item-4746">
          <a href="http://cos.name">
           主页
          </a>
         </li>
         <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-has-children menu-item-8120" id="menu-item-8120">
          <a href="http://cos.name/cn">
           论坛
          </a>
          <ul class="sub-menu">
           <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-8124" id="menu-item-8124">
            <a href="http://cos.name/cn/wp-login.php?action=register">
             论坛注册
            </a>
           </li>
           <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-8125" id="menu-item-8125">
            <a href="http://cos.name/cn/wp-login.php">
             论坛登录
            </a>
           </li>
          </ul>
         </li>
         <li class="menu-item menu-item-type-post_type menu-item-object-page menu-item-8110" id="menu-item-8110">
          <a href="http://cos.name/books/">
           图书资料
          </a>
         </li>
         <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-8112" id="menu-item-8112">
          <a href="http://cos.name/videos">
           视频教程
          </a>
         </li>
         <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-8333" id="menu-item-8333">
          <a href="http://cos.name/salon/">
           统计沙龙
          </a>
         </li>
         <li class="menu-item menu-item-type-post_type menu-item-object-page menu-item-8111" id="menu-item-8111">
          <a href="http://cos.name/chinar/">
           R语言会议
          </a>
         </li>
         <li class="menu-item menu-item-type-post_type menu-item-object-page menu-item-8109" id="menu-item-8109">
          <a href="http://cos.name/training/">
           讲座与培训
          </a>
         </li>
         <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-8115" id="menu-item-8115">
          <a href="http://cos.name/cn/forum/comprehensive/job/">
           招聘信息
          </a>
         </li>
         <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-has-children menu-item-4780" id="menu-item-4780">
          <a href="http://cos.name/about">
           关于我们
          </a>
          <ul class="sub-menu">
           <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-8113" id="menu-item-8113">
            <a href="http://cos.name/2008/11/how-to-work-with-cos/">
             加入我们
            </a>
           </li>
           <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-8119" id="menu-item-8119">
            <a href="http://cos.name/donate/">
             赞助我们
            </a>
           </li>
           <li class="menu-item menu-item-type-custom menu-item-object-custom menu-item-8114" id="menu-item-8114">
            <a href="http://cos.name/cn/forum/26">
             项目合作
            </a>
           </li>
          </ul>
         </li>
        </ul>
       </div>
       <form action="http://cos.name/" class="search-form" method="get" role="search">
        <label>
         <span class="screen-reader-text">
          搜索：
         </span>
        </label>
       </form>
      </nav>
      <!-- #site-navigation -->
     </div>
     <!-- #navbar -->
    </header>
    <!-- #masthead -->
    <div class="site-main" id="main">
     <div class="content-area" id="primary">
      <div class="site-content" id="content" role="main">
       <article class="post-6007 post type-post status-publish format-standard hentry category-optim-sim category-computing category-packages category-software tag-mcmc tag-r-language tag-133 tag-416 tag-244 tag-136 tag-65" id="post-6007">
        <header class="entry-header">
         <h1 class="entry-title">
          MCMC案例学习
         </h1>
         <div class="entry-meta">
          <span class="date">
           <a href="http://cos.name/2012/07/mcmc-case-study/" rel="bookmark" title="链向MCMC案例学习的固定链接">
            <time class="entry-date" datetime="2012-07-24T14:49:32+00:00">
             2012/07/24
            </time>
           </a>
          </span>
          <span class="categories-links">
           <a href="http://cos.name/category/computing/optim-sim/" rel="category tag">
            优化与模拟
           </a>
           、
           <a href="http://cos.name/category/computing/" rel="category tag">
            统计计算
           </a>
           、
           <a href="http://cos.name/category/software/packages/" rel="category tag">
            统计软件
           </a>
           、
           <a href="http://cos.name/category/software/" rel="category tag">
            软件应用
           </a>
          </span>
          <span class="tags-links">
           <a href="http://cos.name/tag/mcmc/" rel="tag">
            MCMC
           </a>
           、
           <a href="http://cos.name/tag/r-language/" rel="tag">
            R语言
           </a>
           、
           <a href="http://cos.name/tag/%e5%8f%82%e6%95%b0%e4%bc%b0%e8%ae%a1/" rel="tag">
            参数估计
           </a>
           、
           <a href="http://cos.name/tag/%e5%b9%bf%e4%b9%89%e7%ba%bf%e6%80%a7%e6%a8%a1%e5%9e%8b/" rel="tag">
            广义线性模型
           </a>
           、
           <a href="http://cos.name/tag/%e6%a8%a1%e6%8b%9f/" rel="tag">
            模拟
           </a>
           、
           <a href="http://cos.name/tag/%e8%ae%a1%e7%ae%97/" rel="tag">
            计算
           </a>
           、
           <a href="http://cos.name/tag/%e8%b4%9d%e5%8f%b6%e6%96%af/" rel="tag">
            贝叶斯
           </a>
          </span>
          <span class="author vcard">
           <a class="url fn n" href="http://cos.name/author/editor/" rel="author" title="查看所有由COS编辑部发布的文章">
            COS编辑部
           </a>
          </span>
         </div>
         <!-- .entry-meta -->
        </header>
        <!-- .entry-header -->
        <div class="entry-content">
         <blockquote>
          <p>
           本文是R中mcmc包的一篇
           <a href="http://www.stat.umn.edu/geyer/mcmc/library/mcmc/doc/demo.pdf" target="_blank">
            帮助文档
           </a>
           ，作者为Charles J.Geyer。经过
           <a href="http://cos.name/2012/06/reproducible-research-with-knitr/" target="_blank">
            knitr
           </a>
           编译后的pdf文档
           <a href="http://cloud.github.com/downloads/cosname/editor/mcmc.pdf" target="_blank">
            可见此处
           </a>
           ，提供中文译稿的作者：
           <br/>
           闫超，天津财经大学统计系2011级研究生，方向：非寿险准备金评估。
           <br/>
           高磊，天津财经大学统计系2011级研究生，方向：非寿险准备金评估。
          </p>
         </blockquote>
         <p>
          这个案例，我们不关心题目的具体意义，重点放在利用贝叶斯的观点来解决问题时，MCMC在后续的计算中所发挥的巨大作用。我们知道，贝叶斯的结果往往是一个后验分布。这个后验分布往往很复杂，我们难以用经典的方法求解其期望与方差等一系列的数据特征，这时MCMC来了，将这一系列问题通过模拟来解决。从这个意义上说，MCMC是一种计算手段。依频率学派看来，题目利用广义线性模型可以解决，在贝叶斯看来同样以解决，但是遇到了一个问题，就是我们得到的非标准后验分布很复杂。我们正是利用MCMC来解决了这个分布的处理问题。本文的重点也在于此。
         </p>
         <p>
          在使用MCMC时作者遵循了这样的思路，首先依照贝叶斯解决问题的套路，构建了非标准后验分布函数。然后初步运行MCMC，确定合适的scale。继而，确定适当的模拟批次和每批长度(以克服模拟取样的相关性)。最后，估计参数并利用delta方法估计标准误。
         </p>
         <h2>
          1. 问题的提出
         </h2>
         <p>
          这是一个关于R软件中
          <strong>
           mcmc
          </strong>
          包的应用案例。问题出自明尼苏达大学统计系博士入学考试试题。这个问题所需要的数据存放在
          <code>
           logit
          </code>
          数据集中。在这个数据集中有五个变量，其中四个自变量
          <code>
           x1、x2、x3、x4
          </code>
          ，一个响应变量
          <code>
           y
          </code>
          。
         </p>
         <p>
          对于这个问题，频率学派的处理方法是利用广义线性模型进行参数估计，下面是相应的R代码以及结果：
         </p>
         <pre><code class="r">library(mcmc)
data(logit)
out &lt;- glm(y ~ x1 + x2 + x3 + x4, data = logit, family = binomial(), x = T)
summary(out)
</code></pre>
         <pre><code>
Call:
glm(formula = y ~ x1 + x2 + x3 + x4, family = binomial(), data = logit, 
    x = T)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.746  -0.691   0.154   0.704   2.194  

Coefficients:
            Estimate Std. Error z value Pr(&gt;|z|)   
(Intercept)    0.633      0.301    2.10   0.0354 * 
x1             0.739      0.362    2.04   0.0410 * 
x2             1.114      0.363    3.07   0.0021 **
x3             0.478      0.354    1.35   0.1766   
x4             0.694      0.399    1.74   0.0817 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 137.628  on 99  degrees of freedom
Residual deviance:  87.668  on 95  degrees of freedom
AIC: 97.67

Number of Fisher Scoring iterations: 6

</code></pre>
         <p>
          但是，使用频率学派的分析方法解决这个问题并不是我们想要的，我们希望使用Bayesian分析方法。对于Bayesian分析而言，我们假定数据模型(即广义线性模型)与频率学派一致。同时假定，五个参数(回归系数)相互独立，并服从均值为0、标准差为2的先验正态分布。
         </p>
         <p>
          定义下面的R函数来计算非标准的对数后验分布概率密度(先验密度与似然函数相乘)。我们为什么要定义成密度函数的对数形式？因为虽然我们是从分布中采样，但是MCMC算法的执行函数
          <code>
           metrop()
          </code>
          需要的一个参数正是这个分布的密度函数的对数形式。
         </p>
         <pre><code class="r">x &lt;- out$x
y &lt;- out$y
lupost &lt;- function(beta, x, y) {
    eta &lt;- as.numeric(x %*% beta)
    logp &lt;- ifelse(eta &lt; 0, eta - log1p(exp(eta)), -log1p(exp(-eta)))
    logq &lt;- ifelse(eta &lt; 0, -log1p(exp(eta)), -eta - log1p(exp(-eta)))
    logl &lt;- sum(logp[y == 1]) + sum(logq[y == 0])
    return(logl - sum(beta^2)/8)
}
</code></pre>
         <p>
          <span id="more-6007">
          </span>
         </p>
         <h2>
          2. 开始MCMC处理
         </h2>
         <p>
          在完成上面数据以及函数的定义(它们都是mcmc算法的输入参数)之后，我们就可以利用
          <strong>
           mcmc
          </strong>
          包中的
          <code>
           metrop()
          </code>
          来产生随机数据，模拟模型的后验分布。也就是说，我们将要从后验分布中进行取样。
         </p>
         <pre><code class="r">set.seed(42)
beta.init &lt;- as.numeric(coefficients(out))
out &lt;- metrop(lupost, beta.init, 1000, x = x, y = y)
names(out)
</code></pre>
         <pre><code> [1] "accept"       "batch"        "initial"      "final"       
 [5] "initial.seed" "final.seed"   "time"         "lud"         
 [9] "nbatch"       "blen"         "nspac"        "scale"       
[13] "debug"       
</code></pre>
         <p>
          此处
          <code>
           metrop()
          </code>
          使用到了如下几种参数：
         </p>
         <ul>
          <li>
           <code>
            lupost
           </code>
           一个函数对象，即后验分布概率密度函数的对数。
          </li>
          <li>
           <code>
            beta.init
           </code>
           用来设定马氏链的初始状态。这里是上面广义线性模型的参数估计结果。
          </li>
          <li>
           <code>
            nbatch = 1000
           </code>
           是采样的样本容量，也就是马氏链的转移次数。
          </li>
          <li>
           <code>
            x、y
           </code>
           ，是传入
           <code>
            lupost
           </code>
           函数中的一些参数。
          </li>
         </ul>
         <p>
          <code>
           metrop()
          </code>
          函数的结果是一个list(列表)。模拟的数据存放在
          <code>
           out$batch
          </code>
          里。刚才的函数运行结果显示接受概率(accept)很低，所以我们调整一下proposal(建议分布)的重要参数
          <code>
           scale
          </code>
          (即随机游走MH算法的方差)来获得一个更合理的接受概率。什么是接受概率？我们用马氏链的方法对状态空间进行采样，那必然我们有一些值是访问不到的，接受概率就是对这种状况进行衡量的。接受概率大，说明访问越充分，接受概率小，访问就受到了一定的限制，访问不是很细致，但样本的自相关性会减弱。因此，接受概率不是越大越好，也不是越小越好。到底多大的接受概率我们认为才是合理的呢？理论显示，在五个需要估计参数的条件下，接受概率达到20%即可。我们开始尝试不同的
          <code>
           scale
          </code>
          值进行测试：
         </p>
         <pre><code class="r">out &lt;- metrop(out, scale = 0.1, x = x, y = y)
out$accept
</code></pre>
         <pre><code>[1] 0.739
</code></pre>
         <pre><code class="r">out &lt;- metrop(out, scale = 0.3, x = x, y = y)
out$accept
</code></pre>
         <pre><code>[1] 0.371
</code></pre>
         <pre><code class="r">out &lt;- metrop(out, scale = 0.5, x = x, y = y)
out$accept
</code></pre>
         <pre><code>[1] 0.148
</code></pre>
         <pre><code class="r">out &lt;- metrop(out, scale = 0.4, x = x, y = y)
out$accept
</code></pre>
         <pre><code>[1] 0.209
</code></pre>
         <pre><code class="r">out &lt;- metrop(out, nbatch = 10000, x = x, y = y)
out$accept
</code></pre>
         <pre><code>[1] 0.2345
</code></pre>
         <p>
          可以看到，每个
          <code>
           metrop()
          </code>
          中第一个参数均为
          <code>
           out
          </code>
          ,它的含义是
          <code>
           metrop()
          </code>
          的许多参数与生成
          <code>
           out
          </code>
          的那一次模拟的参数一致，除了这次又修改的部分外，比如
          <code>
           scale
          </code>
          。这里改变的是
          <code>
           scale
          </code>
          ，这是
          <code>
           Metropolis
          </code>
          随机游走算法中的一个参数值越大，那么马尔可夫链状态转移就越明显(样本自相关性减弱)，但同时接受概率就越低，所以这是一个需要兼顾两者的选择。从运行结果可以看到
          <code>
           scale=0.4
          </code>
          时满足了我们的要求。
         </p>
         <p>
          <!--more-->
         </p>
         <h2>
          3. 诊断，确定批次长度
         </h2>
         <p>
          我们找到了适合接受概率的
          <code>
           scale
          </code>
          值，下面我们按此参数增加模拟次数，然后观察结果。
         </p>
         <pre><code class="r">out &lt;- metrop(out, nbatch = 10000, x = x, y = y)
out$accept
</code></pre>
         <pre><code>[1] 0.2332
</code></pre>
         <pre><code class="r">out$time
</code></pre>
         <pre><code>   user  system elapsed 
   2.40    0.00    2.41 
</code></pre>
         <pre><code class="r">plot(ts(out$batch))
</code></pre>
         <p align="center">
          <img src=""/>
         </p>
         <p>
          解释一下
          <code>
           metrop()
          </code>
          中的参数
          <code>
           nbatch
          </code>
          ，其实与它对应的还有一个参数是
          <code>
           blen
          </code>
          ；
          <code>
           nbatch
          </code>
          是要模拟的批数，
          <code>
           blen
          </code>
          是要模拟的每批的长度，
          <code>
           blen
          </code>
          默认为1。假如
          <code>
           nbatch=100
          </code>
          ，
          <code>
           blen=100
          </code>
          ，那么马氏链共转移100*100=10000次.
         </p>
         <pre><code class="r">acf(out$batch)
</code></pre>
         <p align="center">
          <img src=""/>
         </p>
         <p>
          在这个问题中，我们不研究收敛的问题。从自相关图可以看出，25阶以后的自相关系数可以忽略不计。所以每批次长度25就够了，为了更加保险，我们在下面的模拟中设置每批次长度为100。什么是批次？这是我们为了估计参数而对模拟的结果进行分组时的一个分组长度，之所以要分组，是因为这样才能克服相关性，使批次的均值之间近似独立。这样估计参数才更有效。
         </p>
         <h2>
          4. MCMC 的参数估计值与标准误
         </h2>
         <p>
          首先，运行以下程序：
         </p>
         <pre><code class="r">out &lt;- metrop(out, nbatch = 100, blen = 100, outfun = function(z, ...) c(z, 
    z^2), x = x, y = y)
out$accept
</code></pre>
         <pre><code>[1] 0.2345
</code></pre>
         <pre><code class="r">out$time
</code></pre>
         <pre><code>   user  system elapsed 
   2.46    0.00    2.54 
</code></pre>
         <p>
          <code>
           outfun
          </code>
          函数的作用是方便构建一些统计量。对于这个问题，我们关心的是后验均值和方差。均值的计算相对简单，只需要对涉及的变量进行平均。但是方差的计算有点棘手。我们需要利用等式
          <br/>
          $$
          <br/>
          var(X)=E(X^2)-E(X)^2
          <br/>
          $$
         </p>
         <p>
          将方差表示为可以通过简单的平均进行估计的两部分的方程。因此，我们只需要得到样本的一阶矩和二阶矩。此处，函数
          <code>
           outfun
          </code>
          可针对参数(状态向量)
          <code>
           z
          </code>
          返回
          <code>
           c(z, z
          </code>
          $^2$
          <code>
           )
          </code>
          。
          <code>
           function() c(z, z
          </code>
          $^2$
          <code>
           )
          </code>
          的含义是，每次马氏链转移取样时，得到的一个状态
          <code>
           x
          </code>
          ,把这个状态带入函数中，得到状态本身值，以及它的平方值。这样我们可以求解样本一阶距及二阶矩。
         </p>
         <pre><code class="r">nrow(out$batch)
</code></pre>
         <pre><code>[1] 100
</code></pre>
         <pre><code class="r">out$batch[1, ]
</code></pre>
         <pre><code> [1] 0.6239 0.8217 1.1411 0.4686 0.7494 0.4520 0.7730 1.3696 0.3048 0.6358
</code></pre>
         <p>
          <code>
           out$batch[1, 1]
          </code>
          =
          <code>
           0.6239
          </code>
          是第一批样本的一阶样本矩，
          <code>
           out$batch[1, 6]
          </code>
          =
          <code>
           0.452
          </code>
          是第一批样本的二阶样本矩，注意这里
          <code>
           out$ batch[1, 1]
          </code>
          $^2$与
          <code>
           out$batch[1, 6]
          </code>
          并不相等，因为这里的每批长度(
          <code>
           blen
          </code>
          )是100,只有当
          <code>
           blen=1
          </code>
          时，它们才相等。
         </p>
         <p>
          <code>
           outfun()
          </code>
          函数中参数…是必要的，因为函数同样需要传递其他参数(如这里的
          <code>
           x
          </code>
          和
          <code>
           y
          </code>
          )到
          <code>
           metrop()
          </code>
          。
         </p>
         <h3>
          4.1 简单均值的计算
         </h3>
         <p>
          对每批次的均值再求均值
         </p>
         <pre><code class="r">apply(out$batch, 2, mean)
</code></pre>
         <pre><code> [1] 0.6712 0.7771 1.1814 0.5114 0.7729 0.5465 0.7336 1.5399 0.3878 0.7453
</code></pre>
         <p>
          前五个数就是对后验参数均值的蒙特卡洛估计值，紧随其后的五个数是对后验参数二阶矩的估计值。通过下面的程序，得到参数的方差估计.
         </p>
         <pre><code class="r">foo &lt;- apply(out$batch, 2, mean)
mu &lt;- foo[1:5]
sigmasq &lt;- foo[6:10] - mu^2
mu
</code></pre>
         <pre><code>[1] 0.6712 0.7771 1.1814 0.5114 0.7729
</code></pre>
         <p>
          蒙特卡洛标准误(MCSE)通过批次均值进行计算。这是求均值对简单的方法。
          <br/>
          (注：方差和标准误是两个不同的概念。方差是一个参数估计，而标准误是对参数估计好坏评价的度量。)
         </p>
         <pre><code class="r">mu.muce &lt;- apply(out$batch[, 1:5], 2, sd)/sqrt(out$nbatch)
mu.muce
</code></pre>
         <pre><code>= NA [1] 0.01367 0.01518 0.01785 0.01585 0.01622
</code></pre>
         <p>
          额外因素
          <code>
           sqrt(out$nbatch)
          </code>
          的出现是因为批次均值的方差为$\sigma^{2}/b $，其中b是批次长度，即
          <code>
           out$blen
          </code>
          ；而整体均值的方差为$\sigma^{2}/n$，其中n是迭代总数,即
          <code>
           out$blen
          </code>
          *
          <code>
           out$nbatch
          </code>
          。
         </p>
         <h3>
          4.2 均值的函数
         </h3>
         <p>
          下面我们使用delta method 得到后验方差的MC标准误。$ u_{i} $ 代表某个参数单批次的一阶矩的估计值,$\overline u $代表某个参数所有批次均值的均值，它们都是针对一阶矩而言。对于二阶矩而言样有$ v_{i} $和$ \overline v $ 。令$ u=E(\overline u)$, $ v=E(\overline v) $ 。采用delta方法将非线性函数线性化：
         </p>
         <p>
          $$g(u,v)=v-u^{2}$$
          <br/>
          $$\Delta g(u,v)=\Delta v-2u\Delta u $$
         </p>
         <p>
          也就是说，$ g(\overline u, \overline v)-g(u,v) $ 与 $ (\overline v – v)- 2u(\overline u – u) $具有相同的渐进正态分布。而 $ (\overline v – v)- 2u(\overline u – u) $的方差是$ (v_{i} -v)-2u(u_{i} -u) $方差的1/nbatch倍。这样MCSE可以这样计算：
         </p>
         <p>
          $$\frac{1}{n_{batch}} \displaystyle \sum_{i=1}^{n_{batch}}[(v_{i} – \overline v)-2\overline u
          <br/>
          (u_{i} – \overline u)]^2 $$
         </p>
         <p>
          我们将以上的计算过程用程序实现:
         </p>
         <pre><code class="r">u &lt;- out$batch[, 1:5]
v &lt;- out$batch[, 6:10]
ubar &lt;- apply(u, 2, mean)
vbar &lt;- apply(v, 2, mean)
deltau &lt;- sweep(u, 2, ubar)
deltav &lt;- sweep(v, 2, vbar)
foo &lt;- sweep(deltau, 2, ubar, "*")
sigmasq.mcse &lt;- sqrt(apply((deltav - 2 * foo)^2, 2, mean)/out$nbatch)
sigmasq.mcse
</code></pre>
         <pre><code>[1] 0.004687 0.008586 0.007195 0.007666 0.007714
</code></pre>
         <p>
          这五个值就是后验方差的蒙特卡洛标准误(MCSE)。
         </p>
         <h3>
          4.3 均值函数的函数
         </h3>
         <p>
          如果我们对后验标准差也感兴趣的话，也可以通过delta method计算它的标准误，程序如下
         </p>
         <pre><code class="r">sigma &lt;- sqrt(sigmasq)
sigma.mcse &lt;- sigmasq.mcse/(2 * sigma)
sigma.mcse
</code></pre>
         <pre><code>[1] 0.007565 0.011923 0.009470 0.010789 0.010029
</code></pre>
         <h2>
          5. 最后的运行
         </h2>
         <p>
          问题已经解决。现在唯一需要改进的就是提高结果的精确度。(试题要求“你的马尔科夫链采样器必须足够长以保证参数估计的标准误低于0.01”)。取模拟批次为500，每批长度为400，运行如下：
         </p>
         <pre><code class="r">out &lt;- metrop(out, nbatch = 500, blen = 400, x = x, y = y)
out$accept
</code></pre>
         <pre><code>[1] 0.2352
</code></pre>
         <pre><code class="r">out$time
</code></pre>
         <pre><code>   user  system elapsed 
  50.40    0.01   51.01 
</code></pre>
         <p>
          (显然，由于模拟的次数增大，程序运行时间变长，当然不同运算速度的计算机可能显示结果并不一样。下面进行均值和方差的估计。)
         </p>
         <pre><code class="r">foo &lt;- apply(out$batch, 2, mean)
mu &lt;- foo[1:5]
sigmasq &lt;- foo[6:10] - mu^2
mu
</code></pre>
         <pre><code>[1] 0.6636 0.7961 1.1712 0.5075 0.7241
</code></pre>
         <p>
          然后计算均值估计的标准误
         </p>
         <pre><code class="r">mu.muce &lt;- apply(out$batch[, 1:5], 2, sd)/sqrt(out$nbatch)
mu.muce
</code></pre>
         <pre><code>[1] 0.002972 0.003611 0.003828 0.003604 0.004242
</code></pre>
         <p>
          紧接着计算方差估计的标准误
         </p>
         <pre><code class="r">u &lt;- out$batch[, 1:5]
v &lt;- out$batch[, 6:10]
ubar &lt;- apply(u, 2, mean)
vbar &lt;- apply(v, 2, mean)
deltau &lt;- sweep(u, 2, ubar)
deltav &lt;- sweep(v, 2, vbar)
foo &lt;- sweep(deltau, 2, ubar, "*")
sigmasq.mcse &lt;- sqrt(apply((deltav - 2 * foo)^2, 2, mean)/out$nbatch)
sigmasq.mcse
</code></pre>
         <pre><code>[1] 0.001062 0.001718 0.001538 0.001574 0.002007
</code></pre>
         <p>
          给出标准差的估计以及标准误
         </p>
         <pre><code class="r">sigma &lt;- sqrt(sigmasq)
sigma.mcse &lt;- sigmasq.mcse/(2 * sigma)
sigma.mcse
</code></pre>
         <pre><code>[1] 0.001752 0.002355 0.002116 0.002192 0.002509
</code></pre>
         <p>
          以表格形式展示结果。
         </p>
         <ul>
          <li>
           后验均值估计：
          </li>
         </ul>
         <pre><code class="r">library(xtable)
data1 &lt;- rbind(mu, mu.muce)
colnames(data1) &lt;- c("constant", "x1", "x2", "x3", "x4")
rownames(data1) &lt;- c("estimate", "MCSE")
data1.table &lt;- xtable(data1, digits = 5)
print(data1.table, type = "html")
</code></pre>
         <table border="1">
          <tbody>
           <tr>
            <th>
            </th>
            <th>
             constant
            </th>
            <th>
             x1
            </th>
            <th>
             x2
            </th>
            <th>
             x3
            </th>
            <th>
             x4
            </th>
           </tr>
           <tr>
            <td align="right">
             estimate
            </td>
            <td align="right">
             0.66365
            </td>
            <td align="right">
             0.79608
            </td>
            <td align="right">
             1.17118
            </td>
            <td align="right">
             0.50755
            </td>
            <td align="right">
             0.72414
            </td>
           </tr>
           <tr>
            <td align="right">
             MCSE
            </td>
            <td align="right">
             0.00297
            </td>
            <td align="right">
             0.00361
            </td>
            <td align="right">
             0.00383
            </td>
            <td align="right">
             0.00360
            </td>
            <td align="right">
             0.00424
            </td>
           </tr>
          </tbody>
         </table>
         <ul>
          <li>
           后验方差估计
          </li>
         </ul>
         <pre><code class="r">data2 &lt;- rbind(sigmasq, sigmasq.mcse)
colnames(data2) &lt;- c("constant", "x1", "x2", "x3", "x4")
rownames(data2) &lt;- c("estimate", "MCSE")
data2.table &lt;- xtable(data2, digits = 5)
print(data2.table, type = "html")
</code></pre>
         <table border="1">
          <tbody>
           <tr>
            <th>
            </th>
            <th>
             constant
            </th>
            <th>
             x1
            </th>
            <th>
             x2
            </th>
            <th>
             x3
            </th>
            <th>
             x4
            </th>
           </tr>
           <tr>
            <td align="right">
             estimate
            </td>
            <td align="right">
             0.09195
            </td>
            <td align="right">
             0.13302
            </td>
            <td align="right">
             0.13200
            </td>
            <td align="right">
             0.12890
            </td>
            <td align="right">
             0.15997
            </td>
           </tr>
           <tr>
            <td align="right">
             MCSE
            </td>
            <td align="right">
             0.00106
            </td>
            <td align="right">
             0.00172
            </td>
            <td align="right">
             0.00154
            </td>
            <td align="right">
             0.00157
            </td>
            <td align="right">
             0.00201
            </td>
           </tr>
          </tbody>
         </table>
         <ul>
          <li>
           后验标准差估计
          </li>
         </ul>
         <pre><code class="r">data3 &lt;- rbind(sigma, sigma.mcse)
colnames(data3) &lt;- c("constant", "x1", "x2", "x3", "x4")
rownames(data3) &lt;- c("estimate", "MCSE")
data3.table &lt;- xtable(data3, digits = 5)
print(data3.table, type = "html")
</code></pre>
         <table border="1">
          <tbody>
           <tr>
            <th>
            </th>
            <th>
             constant
            </th>
            <th>
             x1
            </th>
            <th>
             x2
            </th>
            <th>
             x3
            </th>
            <th>
             x4
            </th>
           </tr>
           <tr>
            <td align="right">
             estimate
            </td>
            <td align="right">
             0.30323
            </td>
            <td align="right">
             0.36472
            </td>
            <td align="right">
             0.36332
            </td>
            <td align="right">
             0.35903
            </td>
            <td align="right">
             0.39997
            </td>
           </tr>
           <tr>
            <td align="right">
             MCSE
            </td>
            <td align="right">
             0.00175
            </td>
            <td align="right">
             0.00236
            </td>
            <td align="right">
             0.00212
            </td>
            <td align="right">
             0.00219
            </td>
            <td align="right">
             0.00251
            </td>
           </tr>
          </tbody>
         </table>
         <div class="wumii-hook">
          <br/>
          <br/>
         </div>
        </div>
        <!-- .entry-content -->
        <footer class="entry-meta">
         <div class="author-info">
          <div class="author-avatar">
           <img src="http://sdn.geekzu.org/avatar/2fe058e9e383c85afa949b36e869432f?s=74&amp;d=monsterid&amp;r=g"/>
          </div>
          <!-- .author-avatar -->
          <div class="author-description">
           <h2 class="author-title">
            关于COS编辑部
           </h2>
           <p class="author-bio">
            本账户为COS编辑部公共账户，目前由朱雪宁任主编，王小宁任副主编，编辑有：冯璟烁、吴佳萍、张心雨、施涛、霍志骥、何通、冷静、尤晓斌、肖楠、邱怡轩、高涛、谢益辉等人，主要负责主站文章的规范化编辑以及相关论文、书籍、手册的整理、编纂、出版等工作。
            <a class="author-link" href="http://cos.name/author/editor/" rel="author">
             查看所有由COS编辑部发表的文章
             <span class="meta-nav">
              →
             </span>
            </a>
           </p>
          </div>
          <!-- .author-description -->
         </div>
         <!-- .author-info -->
        </footer>
        <!-- .entry-meta -->
       </article>
       <!-- #post -->
       <nav class="navigation post-navigation" role="navigation">
        <h1 class="screen-reader-text">
         文章导航
        </h1>
        <div class="nav-links">
         <a href="http://cos.name/2012/07/shanghair-july-2012-mindshare/" rel="prev">
          <span class="meta-nav">
           ←
          </span>
          COS数据分析沙龙第二期（上海）
         </a>
         <a href="http://cos.name/2012/08/cos-salon-review-2/" rel="next">
          COS数据分析沙龙第三期（北京）
          <span class="meta-nav">
           →
          </span>
         </a>
        </div>
        <!-- .nav-links -->
       </nav>
       <!-- .navigation -->
       <div class="comments-area" id="comments">
        <h2 class="comments-title">
         《
         <span>
          MCMC案例学习
         </span>
         》有15个想法
        </h2>
        <ol class="comment-list">
         <li class="comment byuser comment-author-dingpeng even thread-even depth-1 parent" id="comment-3238">
          <article class="comment-body" id="div-comment-3238">
           <footer class="comment-meta">
            <div class="comment-author vcard">
             <img src="http://sdn.geekzu.org/avatar/952a7fe51bcad078c1f06cb495ab9851?s=74&amp;d=monsterid&amp;r=g"/>
             <b class="fn">
              丁鹏
             </b>
             <span class="says">
              说道：
             </span>
            </div>
            <!-- .comment-author -->
            <div class="comment-metadata">
             <a href="http://cos.name/2012/07/mcmc-case-study/#comment-3238">
              <time datetime="2012-07-25T08:09:12+00:00">
               2012/07/25 08:09
              </time>
             </a>
            </div>
            <!-- .comment-metadata -->
           </footer>
           <!-- .comment-meta -->
           <div class="comment-content">
            <p>
             拟合bayesian GLM，请用MCMCpack。
            </p>
            <p>
             在现实中，一般不会这样直接的使用metropolis。MCMCpack的算法是正态逼近后验，用来作为metropolis的proposal。
            </p>
            <p>
             如果直接暴力的使用metropolis，大部分时候效果奇差，尤其是X之间的相关性强的时候，mcmc的ACF太大，mixing效果不好。
            </p>
           </div>
           <!-- .comment-content -->
           <div class="reply">
            <a aria-label="回复给丁鹏" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=3238#respond" onclick='return addComment.moveForm( "div-comment-3238", "3238", "respond", "6007" )' rel="nofollow">
             回复
            </a>
           </div>
          </article>
          <!-- .comment-body -->
          <ol class="children">
           <li class="comment odd alt depth-2" id="comment-6457">
            <article class="comment-body" id="div-comment-6457">
             <footer class="comment-meta">
              <div class="comment-author vcard">
               <img src="http://sdn.geekzu.org/avatar/bf9c96b3196a7123890955ed98325c63?s=74&amp;d=monsterid&amp;r=g"/>
               <b class="fn">
                akward
               </b>
               <span class="says">
                说道：
               </span>
              </div>
              <!-- .comment-author -->
              <div class="comment-metadata">
               <a href="http://cos.name/2012/07/mcmc-case-study/#comment-6457">
                <time datetime="2014-12-17T21:41:02+00:00">
                 2014/12/17 21:41
                </time>
               </a>
              </div>
              <!-- .comment-metadata -->
             </footer>
             <!-- .comment-meta -->
             <div class="comment-content">
              <p>
               请问，我的代码：
               <br/>
               library(MCMCpack)
               <br/>
               out1 &lt;- MCMCprobit(y~et+TFP+Size+Wage+Age+Mid,data=export_data,burnin = 5000, mcmc = 10000,
               <br/>
               b0 = c(0,0,0,0,0,0), B0 = c(100,100,100,100,100,100), marginal.likelihood="Chib95")
               <br/>
               plot(out2)
               <br/>
               summary(out2)
              </p>
              <p>
               参数估计结果里有个sigma是什么？
               <br/>
               这个是用来衡量什么的？
               <br/>
               模型估计的好坏看什么？有拟合优度指标没？
              </p>
             </div>
             <!-- .comment-content -->
             <div class="reply">
              <a aria-label="回复给akward" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=6457#respond" onclick='return addComment.moveForm( "div-comment-6457", "6457", "respond", "6007" )' rel="nofollow">
               回复
              </a>
             </div>
            </article>
            <!-- .comment-body -->
           </li>
           <!-- #comment-## -->
           <li class="comment even depth-2" id="comment-7130">
            <article class="comment-body" id="div-comment-7130">
             <footer class="comment-meta">
              <div class="comment-author vcard">
               <img src="http://sdn.geekzu.org/avatar/3c8033452a17b9d138179e4a0c56cb95?s=74&amp;d=monsterid&amp;r=g"/>
               <b class="fn">
                陈宇
               </b>
               <span class="says">
                说道：
               </span>
              </div>
              <!-- .comment-author -->
              <div class="comment-metadata">
               <a href="http://cos.name/2012/07/mcmc-case-study/#comment-7130">
                <time datetime="2016-04-02T16:54:31+00:00">
                 2016/04/02 16:54
                </time>
               </a>
              </div>
              <!-- .comment-metadata -->
             </footer>
             <!-- .comment-meta -->
             <div class="comment-content">
              <p>
               大神，如果我只是相对均匀分布 和 正态分布进行采样，能用mcmc方法嘛
              </p>
             </div>
             <!-- .comment-content -->
             <div class="reply">
              <a aria-label="回复给陈宇" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=7130#respond" onclick='return addComment.moveForm( "div-comment-7130", "7130", "respond", "6007" )' rel="nofollow">
               回复
              </a>
             </div>
            </article>
            <!-- .comment-body -->
           </li>
           <!-- #comment-## -->
          </ol>
          <!-- .children -->
         </li>
         <!-- #comment-## -->
         <li class="comment odd alt thread-odd thread-alt depth-1 parent" id="comment-3239">
          <article class="comment-body" id="div-comment-3239">
           <footer class="comment-meta">
            <div class="comment-author vcard">
             <img src="http://sdn.geekzu.org/avatar/5302fa01fd9ecc79a9b95fa6e71a4b19?s=74&amp;d=monsterid&amp;r=g"/>
             <b class="fn">
              gaolei
             </b>
             <span class="says">
              说道：
             </span>
            </div>
            <!-- .comment-author -->
            <div class="comment-metadata">
             <a href="http://cos.name/2012/07/mcmc-case-study/#comment-3239">
              <time datetime="2012-07-25T21:25:40+00:00">
               2012/07/25 21:25
              </time>
             </a>
            </div>
            <!-- .comment-metadata -->
           </footer>
           <!-- .comment-meta -->
           <div class="comment-content">
            <p>
             的确，处理实际问题直接使用这个函数确实不妥。其实，从贝叶斯角度解决广义线性模型，也有这样的包和函数，直接来做，比如，arm 包中的bayesglm(),与glm的使用差不多，不过其实质还是mcmc。处理各种模型，使用贝叶斯方法很普遍，参见（CRAN Task View: Bayesian Inference）。这里只是向初学者展示一下mcmc的使用过程。
            </p>
           </div>
           <!-- .comment-content -->
           <div class="reply">
            <a aria-label="回复给gaolei" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=3239#respond" onclick='return addComment.moveForm( "div-comment-3239", "3239", "respond", "6007" )' rel="nofollow">
             回复
            </a>
           </div>
          </article>
          <!-- .comment-body -->
          <ol class="children">
           <li class="comment even depth-2" id="comment-7131">
            <article class="comment-body" id="div-comment-7131">
             <footer class="comment-meta">
              <div class="comment-author vcard">
               <img src="http://sdn.geekzu.org/avatar/3c8033452a17b9d138179e4a0c56cb95?s=74&amp;d=monsterid&amp;r=g"/>
               <b class="fn">
                陈宇
               </b>
               <span class="says">
                说道：
               </span>
              </div>
              <!-- .comment-author -->
              <div class="comment-metadata">
               <a href="http://cos.name/2012/07/mcmc-case-study/#comment-7131">
                <time datetime="2016-04-02T16:56:07+00:00">
                 2016/04/02 16:56
                </time>
               </a>
              </div>
              <!-- .comment-metadata -->
             </footer>
             <!-- .comment-meta -->
             <div class="comment-content">
              <p>
               大神，如果我只是相对 均匀分布 和 正态分布进行采样，您觉得可行嘛？那建议分布就是正态分布本身嘛？
              </p>
             </div>
             <!-- .comment-content -->
             <div class="reply">
              <a aria-label="回复给陈宇" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=7131#respond" onclick='return addComment.moveForm( "div-comment-7131", "7131", "respond", "6007" )' rel="nofollow">
               回复
              </a>
             </div>
            </article>
            <!-- .comment-body -->
           </li>
           <!-- #comment-## -->
          </ol>
          <!-- .children -->
         </li>
         <!-- #comment-## -->
         <li class="comment odd alt thread-even depth-1" id="comment-3240">
          <article class="comment-body" id="div-comment-3240">
           <footer class="comment-meta">
            <div class="comment-author vcard">
             <img src="http://sdn.geekzu.org/avatar/283d380eb3180db41eabfc919a3655f8?s=74&amp;d=monsterid&amp;r=g"/>
             <b class="fn">
              Lu.M
             </b>
             <span class="says">
              说道：
             </span>
            </div>
            <!-- .comment-author -->
            <div class="comment-metadata">
             <a href="http://cos.name/2012/07/mcmc-case-study/#comment-3240">
              <time datetime="2012-07-25T21:32:26+00:00">
               2012/07/25 21:32
              </time>
             </a>
            </div>
            <!-- .comment-metadata -->
           </footer>
           <!-- .comment-meta -->
           <div class="comment-content">
            <p>
             太好了，虚心学习。二位继续努力。
            </p>
           </div>
           <!-- .comment-content -->
           <div class="reply">
            <a aria-label="回复给Lu.M" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=3240#respond" onclick='return addComment.moveForm( "div-comment-3240", "3240", "respond", "6007" )' rel="nofollow">
             回复
            </a>
           </div>
          </article>
          <!-- .comment-body -->
         </li>
         <!-- #comment-## -->
         <li class="pingback even thread-odd thread-alt depth-1" id="comment-3244">
          <div class="comment-body">
           Pingback：
           <a class="url" href="http://r-ke.info/2012/07/27/mcmc-case-study.html" rel="external nofollow">
            cos: MCMC案例学习 | R客 - 新鲜的中文R语言资讯
           </a>
          </div>
         </li>
         <!-- #comment-## -->
         <li class="comment odd alt thread-even depth-1 parent" id="comment-3246">
          <article class="comment-body" id="div-comment-3246">
           <footer class="comment-meta">
            <div class="comment-author vcard">
             <img src="http://sdn.geekzu.org/avatar/4ff82f7841e42f0a7b1e398e4dc284eb?s=74&amp;d=monsterid&amp;r=g"/>
             <b class="fn">
              <a class="url" href="http://www.minidata.org/" rel="external nofollow">
               白日梦
              </a>
             </b>
             <span class="says">
              说道：
             </span>
            </div>
            <!-- .comment-author -->
            <div class="comment-metadata">
             <a href="http://cos.name/2012/07/mcmc-case-study/#comment-3246">
              <time datetime="2012-07-28T09:29:47+00:00">
               2012/07/28 09:29
              </time>
             </a>
            </div>
            <!-- .comment-metadata -->
           </footer>
           <!-- .comment-meta -->
           <div class="comment-content">
            <p>
             能不能在开篇的时候多讲讲建模原理和步骤, 一上来就是R代码, R基础不好的人看的不是很明白啊.
            </p>
           </div>
           <!-- .comment-content -->
           <div class="reply">
            <a aria-label="回复给白日梦" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=3246#respond" onclick='return addComment.moveForm( "div-comment-3246", "3246", "respond", "6007" )' rel="nofollow">
             回复
            </a>
           </div>
          </article>
          <!-- .comment-body -->
          <ol class="children">
           <li class="comment byuser comment-author-yihui even depth-2" id="comment-3253">
            <article class="comment-body" id="div-comment-3253">
             <footer class="comment-meta">
              <div class="comment-author vcard">
               <img src="http://sdn.geekzu.org/avatar/1022d8e6ebc94e8f6bca9a86cebe312a?s=74&amp;d=monsterid&amp;r=g"/>
               <b class="fn">
                <a class="url" href="http://yihui.name" rel="external nofollow">
                 谢益辉
                </a>
               </b>
               <span class="says">
                说道：
               </span>
              </div>
              <!-- .comment-author -->
              <div class="comment-metadata">
               <a href="http://cos.name/2012/07/mcmc-case-study/#comment-3253">
                <time datetime="2012-07-28T12:11:50+00:00">
                 2012/07/28 12:11
                </time>
               </a>
              </div>
              <!-- .comment-metadata -->
             </footer>
             <!-- .comment-meta -->
             <div class="comment-content">
              <p>
               这本来就是R包的文档，不可避免要上R。一篇文章不可能面面俱到，更多的功课还得你自己去做，况且无代码无真相啊。
              </p>
             </div>
             <!-- .comment-content -->
             <div class="reply">
              <a aria-label="回复给谢益辉" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=3253#respond" onclick='return addComment.moveForm( "div-comment-3253", "3253", "respond", "6007" )' rel="nofollow">
               回复
              </a>
             </div>
            </article>
            <!-- .comment-body -->
           </li>
           <!-- #comment-## -->
           <li class="comment odd alt depth-2" id="comment-3260">
            <article class="comment-body" id="div-comment-3260">
             <footer class="comment-meta">
              <div class="comment-author vcard">
               <img src="http://sdn.geekzu.org/avatar/5302fa01fd9ecc79a9b95fa6e71a4b19?s=74&amp;d=monsterid&amp;r=g"/>
               <b class="fn">
                gaolei
               </b>
               <span class="says">
                说道：
               </span>
              </div>
              <!-- .comment-author -->
              <div class="comment-metadata">
               <a href="http://cos.name/2012/07/mcmc-case-study/#comment-3260">
                <time datetime="2012-07-31T15:53:25+00:00">
                 2012/07/31 15:53
                </time>
               </a>
              </div>
              <!-- .comment-metadata -->
             </footer>
             <!-- .comment-meta -->
             <div class="comment-content">
              <p>
               一块是贝叶斯推断，一块是MCMC算法，看一下这方面的文档就可以了
              </p>
             </div>
             <!-- .comment-content -->
             <div class="reply">
              <a aria-label="回复给gaolei" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=3260#respond" onclick='return addComment.moveForm( "div-comment-3260", "3260", "respond", "6007" )' rel="nofollow">
               回复
              </a>
             </div>
            </article>
            <!-- .comment-body -->
           </li>
           <!-- #comment-## -->
          </ol>
          <!-- .children -->
         </li>
         <!-- #comment-## -->
         <li class="comment even thread-odd thread-alt depth-1" id="comment-4152">
          <article class="comment-body" id="div-comment-4152">
           <footer class="comment-meta">
            <div class="comment-author vcard">
             <img src="http://sdn.geekzu.org/avatar/1ff898f0a9eea1f9ad819f49448ee277?s=74&amp;d=monsterid&amp;r=g"/>
             <b class="fn">
              zhudan
             </b>
             <span class="says">
              说道：
             </span>
            </div>
            <!-- .comment-author -->
            <div class="comment-metadata">
             <a href="http://cos.name/2012/07/mcmc-case-study/#comment-4152">
              <time datetime="2013-05-02T21:07:50+00:00">
               2013/05/02 21:07
              </time>
             </a>
            </div>
            <!-- .comment-metadata -->
           </footer>
           <!-- .comment-meta -->
           <div class="comment-content">
            <p>
             虚心学习！
            </p>
           </div>
           <!-- .comment-content -->
           <div class="reply">
            <a aria-label="回复给zhudan" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=4152#respond" onclick='return addComment.moveForm( "div-comment-4152", "4152", "respond", "6007" )' rel="nofollow">
             回复
            </a>
           </div>
          </article>
          <!-- .comment-body -->
         </li>
         <!-- #comment-## -->
         <li class="comment odd alt thread-even depth-1" id="comment-4179">
          <article class="comment-body" id="div-comment-4179">
           <footer class="comment-meta">
            <div class="comment-author vcard">
             <img src="http://sdn.geekzu.org/avatar/d64c4947e5e63dab06b610ba7326b871?s=74&amp;d=monsterid&amp;r=g"/>
             <b class="fn">
              xwy
             </b>
             <span class="says">
              说道：
             </span>
            </div>
            <!-- .comment-author -->
            <div class="comment-metadata">
             <a href="http://cos.name/2012/07/mcmc-case-study/#comment-4179">
              <time datetime="2013-05-09T23:39:38+00:00">
               2013/05/09 23:39
              </time>
             </a>
            </div>
            <!-- .comment-metadata -->
           </footer>
           <!-- .comment-meta -->
           <div class="comment-content">
            <p>
             这应该是对mcmc package demo的翻译，值得学习
            </p>
           </div>
           <!-- .comment-content -->
           <div class="reply">
            <a aria-label="回复给xwy" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=4179#respond" onclick='return addComment.moveForm( "div-comment-4179", "4179", "respond", "6007" )' rel="nofollow">
             回复
            </a>
           </div>
          </article>
          <!-- .comment-body -->
         </li>
         <!-- #comment-## -->
         <li class="comment even thread-odd thread-alt depth-1 parent" id="comment-4949">
          <article class="comment-body" id="div-comment-4949">
           <footer class="comment-meta">
            <div class="comment-author vcard">
             <img src="http://sdn.geekzu.org/avatar/5241bb4d973a0eaade0713f5f0267ac2?s=74&amp;d=monsterid&amp;r=g"/>
             <b class="fn">
              十八学士
             </b>
             <span class="says">
              说道：
             </span>
            </div>
            <!-- .comment-author -->
            <div class="comment-metadata">
             <a href="http://cos.name/2012/07/mcmc-case-study/#comment-4949">
              <time datetime="2013-10-08T21:20:32+00:00">
               2013/10/08 21:20
              </time>
             </a>
            </div>
            <!-- .comment-metadata -->
           </footer>
           <!-- .comment-meta -->
           <div class="comment-content">
            <p>
             如果其中有个参数有限制，比如非负，可以让这个参数先验分布为一个对数正态分布吗？
            </p>
           </div>
           <!-- .comment-content -->
           <div class="reply">
            <a aria-label="回复给十八学士" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=4949#respond" onclick='return addComment.moveForm( "div-comment-4949", "4949", "respond", "6007" )' rel="nofollow">
             回复
            </a>
           </div>
          </article>
          <!-- .comment-body -->
          <ol class="children">
           <li class="comment odd alt depth-2" id="comment-5181">
            <article class="comment-body" id="div-comment-5181">
             <footer class="comment-meta">
              <div class="comment-author vcard">
               <img src="http://sdn.geekzu.org/avatar/?s=74&amp;d=monsterid&amp;r=g"/>
               <b class="fn">
                <a class="url" href="http://weibo.com/1374894974" rel="external nofollow">
                 lei_gao
                </a>
               </b>
               <span class="says">
                说道：
               </span>
              </div>
              <!-- .comment-author -->
              <div class="comment-metadata">
               <a href="http://cos.name/2012/07/mcmc-case-study/#comment-5181">
                <time datetime="2013-12-15T21:33:13+00:00">
                 2013/12/15 21:33
                </time>
               </a>
              </div>
              <!-- .comment-metadata -->
             </footer>
             <!-- .comment-meta -->
             <div class="comment-content">
              <p>
               没错，用对数正态分布作为其先验分布，能把参数限制在正数上。
              </p>
             </div>
             <!-- .comment-content -->
             <div class="reply">
              <a aria-label="回复给lei_gao" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=5181#respond" onclick='return addComment.moveForm( "div-comment-5181", "5181", "respond", "6007" )' rel="nofollow">
               回复
              </a>
             </div>
            </article>
            <!-- .comment-body -->
           </li>
           <!-- #comment-## -->
          </ol>
          <!-- .children -->
         </li>
         <!-- #comment-## -->
         <li class="comment even thread-even depth-1" id="comment-6505">
          <article class="comment-body" id="div-comment-6505">
           <footer class="comment-meta">
            <div class="comment-author vcard">
             <img src="http://sdn.geekzu.org/avatar/?s=74&amp;d=monsterid&amp;r=g"/>
             <b class="fn">
              <a class="url" href="http://t.qq.com/TEUFELCAT" rel="external nofollow">
               第九街TEUFEL
              </a>
             </b>
             <span class="says">
              说道：
             </span>
            </div>
            <!-- .comment-author -->
            <div class="comment-metadata">
             <a href="http://cos.name/2012/07/mcmc-case-study/#comment-6505">
              <time datetime="2015-01-14T01:22:12+00:00">
               2015/01/14 01:22
              </time>
             </a>
            </div>
            <!-- .comment-metadata -->
           </footer>
           <!-- .comment-meta -->
           <div class="comment-content">
            <p>
             请问，有没有人能告诉我，metrop（lupost, beta.init, 1000, x = x, y = y)，这个写法其中x=x，y=y是定义lupost的参数x，y为样本x，y么？假设lupost的样本变量只有y，样本信息也只有一个变量y，此处应该怎样写呢？
            </p>
           </div>
           <!-- .comment-content -->
           <div class="reply">
            <a aria-label="回复给第九街TEUFEL" class="comment-reply-link" href="http://cos.name/2012/07/mcmc-case-study/?replytocom=6505#respond" onclick='return addComment.moveForm( "div-comment-6505", "6505", "respond", "6007" )' rel="nofollow">
             回复
            </a>
           </div>
          </article>
          <!-- .comment-body -->
         </li>
         <!-- #comment-## -->
        </ol>
        <!-- .comment-list -->
        <div class="comment-respond" id="respond">
         <h3 class="comment-reply-title" id="reply-title">
          发表评论
          <small>
           <a href="http://cos.name/2012/07/mcmc-case-study/#respond" id="cancel-comment-reply-link" rel="nofollow" style="display:none;">
            取消回复
           </a>
          </small>
         </h3>
         <form action="http://cos.name/wp-comments-post.php" class="comment-form" id="commentform" method="post" novalidate="">
          <p class="comment-notes">
           <span id="email-notes">
            电子邮件地址不会被公开。
           </span>
           必填项已用
           <span class="required">
            *
           </span>
           标注
          </p>
          <p class="comment-form-comment">
           <label for="comment">
            评论
           </label>
           <textarea aria-required="true" cols="45" id="comment" maxlength="65525" name="comment" required="required" rows="8">
           </textarea>
          </p>
          <p class="comment-form-author">
           <label for="author">
            姓名
            <span class="required">
             *
            </span>
           </label>
          </p>
          <p class="comment-form-email">
           <label for="email">
            电子邮件
            <span class="required">
             *
            </span>
           </label>
          </p>
          <p class="comment-form-url">
           <label for="url">
            站点
           </label>
          </p>
          <p class="form-submit">
          </p>
          <p style="display: none;">
          </p>
          <p style="display: none;">
          </p>
         </form>
        </div>
        <!-- #respond -->
       </div>
       <!-- #comments -->
      </div>
      <!-- #content -->
     </div>
     <!-- #primary -->
     <div class="sidebar-container" id="tertiary" role="complementary">
      <div class="sidebar-inner">
       <div class="widget-area">
        <aside class="widget widget_text" id="text-3">
         <h3 class="widget-title">
          关注统计之都
         </h3>
         <div class="textwidget">
          <ul>
           <li>
            新浪微博
            <a href="http://weibo.com/cosname">
             @统计之都
            </a>
           </li>
           <li>
            人人网
            <a href="http://renren.com/cosname">
             @统计之都
            </a>
           </li>
           <li>
            Twitter
            <a href="http://twitter.com/cos_name">
             @cos_name
            </a>
           </li>
          </ul>
         </div>
        </aside>
        <aside class="widget widget_text" id="text-6">
         <h3 class="widget-title">
          微信公众平台
         </h3>
         <div class="textwidget">
          <img src="http://cos.name/wp-content/uploads/2013/04/qrcode-8cm.jpg"/>
          <p style="font-size:12px;margin-left:15px">
           微信号 CapStat
          </p>
          <p>
           我们将第一时间向您推送主站和论坛的精彩内容，以及统计之都的线下活动、竞赛、培训和会议信息。
          </p>
         </div>
        </aside>
        <aside class="widget widget_text" id="text-8">
         <h3 class="widget-title">
          站内导航
         </h3>
         <div class="textwidget">
          <ul>
           <li class="page_item">
            <a href="http://cos.name/cn/">
             中文论坛
            </a>
           </li>
           <li class="page_item">
            <a href="http://cos.name/videos/">
             视频
            </a>
           </li>
           <li class="page_item">
            <a href="http://cos.name/salon/">
             统计沙龙
            </a>
           </li>
           <li class="page_item">
            <a href="http://cos.name/books/">
             图书出版
            </a>
           </li>
           <li class="page_item">
            <a href="http://cos.name/training/">
             教育培训
            </a>
           </li>
           <li class="page_item">
            <a href="http://cos.name/about/">
             关于我们
            </a>
           </li>
           <li class="page_item">
            <a href="http://cos.name/donate/">
             捐赠
            </a>
           </li>
           <li class="page_item">
            <a href="http://cos.name/chinar/">
             R语言会议
            </a>
           </li>
          </ul>
         </div>
        </aside>
        <aside class="widget widget_links" id="linkcat-2">
         <h3 class="widget-title">
          友情链接
         </h3>
         <ul class="xoxo blogroll">
          <li>
           <a href="http://stat.ruc.edu.cn" target="_blank" title="中国人民大学统计学院网站">
            中国人民大学统计学院
           </a>
          </li>
          <li>
           <a href="http://rucdmc.net">
            中国人民大学数据挖掘中心
           </a>
          </li>
          <li>
           <a href="http://birc.gsm.pku.edu.cn/" target="_blank">
            北京大学商务智能研究中心
           </a>
          </li>
          <li>
           <a href="http://sam.cufe.edu.cn/" target="_blank" title="中央财经大学统计与数学学院网站">
            中央财经大学统计与数学学院
           </a>
          </li>
          <li>
           <a href="http://tjx.cueb.edu.cn/" target="_blank" title="首都经济贸易大学统计学院网站">
            首经贸统计学院
           </a>
          </li>
          <li>
           <a href="http://www.shookr.com/">
            数客网大数据社区
           </a>
          </li>
          <li>
           <a href="http://www.xueqing.tv/" target="_blank" title="数据科学在线学习平台">
            雪晴数据网
           </a>
          </li>
          <li>
           <a href="http://iera.name/" target="_blank" title="IERA是一个旨在普及、传播和增进工业工程知识的非营利性网站">
            IERA（直通IE）
           </a>
          </li>
         </ul>
        </aside>
        <aside class="widget widget_categories" id="categories-2">
         <h3 class="widget-title">
          全部分类
         </h3>
         <label class="screen-reader-text" for="cat">
          全部分类
         </label>
         <select class="postform" id="cat" name="cat">
          <option value="-1">
           选择分类目录
          </option>
          <option class="level-0" value="925">
           cos访谈  (4)
          </option>
          <option class="level-0" value="659">
           图书出版  (5)
          </option>
          <option class="level-0" value="379">
           数学方法  (14)
          </option>
          <option class="level-1" value="381">
           分析与代数  (1)
          </option>
          <option class="level-1" value="380">
           概率论  (9)
          </option>
          <option class="level-1" value="382">
           随机过程  (5)
          </option>
          <option class="level-0" value="210">
           数据分析  (81)
          </option>
          <option class="level-1" value="203">
           多元统计  (3)
          </option>
          <option class="level-1" value="42">
           数据挖掘与机器学习  (42)
          </option>
          <option class="level-1" value="36">
           生物与医学统计  (17)
          </option>
          <option class="level-1" value="35">
           计量经济学  (4)
          </option>
          <option class="level-1" value="296">
           金融统计  (3)
          </option>
          <option class="level-1" value="303">
           风险精算  (7)
          </option>
          <option class="level-0" value="177">
           模型专题  (15)
          </option>
          <option class="level-1" value="38">
           回归分析  (10)
          </option>
          <option class="level-1" value="41">
           时间序列  (2)
          </option>
          <option class="level-0" value="784">
           每周精选  (24)
          </option>
          <option class="level-1" value="183">
           可视化  (9)
          </option>
          <option class="level-0" value="967">
           沙龙纪要  (3)
          </option>
          <option class="level-0" value="18">
           经典理论  (46)
          </option>
          <option class="level-1" value="37">
           抽样调查  (3)
          </option>
          <option class="level-1" value="4">
           统计推断  (26)
          </option>
          <option class="level-1" value="236">
           试验设计  (7)
          </option>
          <option class="level-1" value="39">
           非参数统计  (3)
          </option>
          <option class="level-0" value="1">
           统计之都  (279)
          </option>
          <option class="level-1" value="884">
           中国R会议  (2)
          </option>
          <option class="level-1" value="885">
           中国R语言会议  (5)
          </option>
          <option class="level-1" value="446">
           出国留学  (3)
          </option>
          <option class="level-1" value="179">
           推荐文章  (90)
          </option>
          <option class="level-1" value="3">
           新闻通知  (75)
          </option>
          <option class="level-1" value="263">
           统计刊物  (10)
          </option>
          <option class="level-1" value="174">
           网站导读  (40)
          </option>
          <option class="level-1" value="204">
           职业事业  (51)
          </option>
          <option class="level-1" value="213">
           高校课堂  (9)
          </option>
          <option class="level-0" value="178">
           统计计算  (28)
          </option>
          <option class="level-1" value="40">
           优化与模拟  (15)
          </option>
          <option class="level-1" value="43">
           贝叶斯方法  (6)
          </option>
          <option class="level-0" value="378">
           软件应用  (116)
          </option>
          <option class="level-1" value="44">
           统计图形  (36)
          </option>
          <option class="level-1" value="110">
           统计软件  (83)
          </option>
         </select>
        </aside>
        <aside class="widget widget_recent_comments" id="recent-comments-2">
         <h3 class="widget-title">
          最新评论
         </h3>
         <ul id="recentcomments">
          <li class="recentcomments">
           <span class="comment-author-link">
            fineboom
           </span>
           发表在《
           <a href="http://cos.name/2016/06/use-shiny-fleetly-set-up-visual-prototype-system/#comment-7317">
            利用shiny包快速搭建可视化原型系统
           </a>
           》
          </li>
          <li class="recentcomments">
           <span class="comment-author-link">
            胡家新
           </span>
           发表在《
           <a href="http://cos.name/2016/06/r%e8%af%ad%e5%8d%83%e5%af%bb%e7%ac%ac%e4%b8%89%e6%9c%9f%ef%bc%9a%e5%bc%a0%e6%97%a0%e5%bf%8c%e7%a9%b6%e7%ab%9f%e7%88%b1%e8%b0%81%ef%bc%9f/#comment-7316">
            R语千寻第三期：张无忌究竟爱谁？
           </a>
           》
          </li>
          <li class="recentcomments">
           <span class="comment-author-link">
            <a class="url" href="http://www.zijiacha.com/category.php?id=6" rel="external nofollow">
             南糯山普洱茶
            </a>
           </span>
           发表在《
           <a href="http://cos.name/2016/06/r%e8%af%ad%e5%8d%83%e5%af%bb%e7%ac%ac%e4%b8%89%e6%9c%9f%ef%bc%9a%e5%bc%a0%e6%97%a0%e5%bf%8c%e7%a9%b6%e7%ab%9f%e7%88%b1%e8%b0%81%ef%bc%9f/#comment-7315">
            R语千寻第三期：张无忌究竟爱谁？
           </a>
           》
          </li>
          <li class="recentcomments">
           <span class="comment-author-link">
            J
           </span>
           发表在《
           <a href="http://cos.name/2016/05/value-of-the-reputation-from-the-data/#comment-7314">
            数据告诉你：高信誉的卖家应该收高价，还是收低价？
           </a>
           》
          </li>
          <li class="recentcomments">
           <span class="comment-author-link">
            <a class="url" href="http://gg" rel="external nofollow">
             Hilda
            </a>
           </span>
           发表在《
           <a href="http://cos.name/2013/01/drawing-map-in-r-era/#comment-7311">
            R时代，你要怎样画地图？
           </a>
           》
          </li>
         </ul>
        </aside>
        <aside class="widget widget_rss" id="rss-282869971">
         <h3 class="widget-title">
          <a class="rsswidget" href="http://cos.name/cn/topics/feed/">
           <img src="http://cos.name/wp-includes/images/rss.png"/>
          </a>
          <a class="rsswidget" href="http://cos.name/cn/topics/feed/">
           中文论坛新帖
          </a>
         </h3>
         <ul>
          <li>
           <a class="rsswidget" href="http://cos.name/cn/topic/417372/">
            处理时间数据和产生时间序列的问题
           </a>
          </li>
          <li>
           <a class="rsswidget" href="http://cos.name/cn/topic/6790/">
            《统计陷阱》下载 （How to lie with statistics）
           </a>
          </li>
          <li>
           <a class="rsswidget" href="http://cos.name/cn/topic/818/">
            统计学的世界（第五版）
           </a>
          </li>
          <li>
           <a class="rsswidget" href="http://cos.name/cn/topic/16574/">
            class(x) 返回值值是AsIs,AsIs代表什么，有什么用处？
           </a>
          </li>
          <li>
           <a class="rsswidget" href="http://cos.name/cn/topic/417366/">
            如何用R绘制一个分类算法的决策规则
           </a>
          </li>
         </ul>
        </aside>
        <aside class="widget widget_text" id="text-2">
         <h3 class="widget-title">
          登录/RSS
         </h3>
         <div class="textwidget">
          <ul>
           <li>
            <a href="http://cos.name/wp-admin/">
             登录
            </a>
           </li>
           <li>
            <a href="http://cos.name/feed/" title="使用 RSS 2.0 同步站点内容">
             文章
             <abbr title="Really Simple Syndication">
              RSS
             </abbr>
            </a>
           </li>
           <li>
            <a href="http://cos.name/comments/feed/" title="RSS 上的最近评论">
             评论
             <abbr title="Really Simple Syndication">
              RSS
             </abbr>
            </a>
           </li>
          </ul>
         </div>
        </aside>
       </div>
       <!-- .widget-area -->
      </div>
      <!-- .sidebar-inner -->
     </div>
     <!-- #tertiary -->
    </div>
    <!-- #main -->
    <footer class="site-footer" id="colophon" role="contentinfo">
     <div class="site-info">
      版权所有 © 2014 统计之都 | 由
      <a href="http://wordpress.org/">
       WordPress
      </a>
      构建 | 主题修改自
      <a href="http://wordpress.org/themes/twentythirteen">
       Twenty Thirteen
      </a>
     </div>
     <!-- .site-info -->
    </footer>
    <!-- #colophon -->
   </div>
   <!-- #page -->
   <p style="margin:0;padding:0;height:1px;overflow:hidden;">
    <a href="http://www.wumii.com/widget/relatedItems" style="border:0;">
     <img src="http://static.wumii.cn/images/pixel.png"/>
    </a>
   </p>
  </div>
 </body>
</html>